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Abstract. We address phase-shift estimation by means of squeezed vacuum 
probe and homodyne detection. We analyze Bayesian estimator, which is known 
to asymptotically saturate the classical Cramer-Rao bound to the variance, and 
discuss convergence looking at the a posteriori distribution as the number of 
measurements increases. We also suggest two feasible adaptive methods, acting on 
the squeezing parameter and/or the homodyne local oscillator phase, which allow 
to optimize homodyne detection and approach the ultimate bound to precision 
imposed by the quantum Cramer- Rao theorem. The performances of our two-step 
methods are investigated by means of Monte Carlo simulated experiments with a 
small number of homodyne data, thus giving a quantitative meaning to the notion 
of asymptotic optimality. 



Bayesian estimation in homodyne interferometry 



2 



1. Introduction 

Quantum phase measurements cannot be described by means of a proper observable 
and different operational approaches have been introduced over the years [H [21 [3l [H 
[51 [3] . On the other hand, from a practical point of view, phase detection of quantum 
fields is generally associated with interferometric devices, i.e., detection schemes aimed 
at the estimation of phase by measuring field- or intensity-based quantities with phase- 
dependent statistics [71 [HI [SI [101 [HI HH H3] • The art of interferomety, in turn, consists 
in answering to two question: a) How can the unknown phase be effectively retrieved 
from the data sample? and b) Which is the resulting precision? The first point 
amounts to the choice of an estimator, i.e., a function from the data sample to the set of 
possible values of the phase-shift. Among possible estimators Bayes p~4] and maximum 
likelihood ones [HI [16] play a special role due to their asympotic (i.e., for large 
number of measurements) properties. The second point may be properly addressed 
in the framework of quantum estimation theory, which addresses the inference of a 
physical quantity which is not directly accessible by means of the measurement of 
a different observable, or a set of observables, somehow related to the quantity of 
interest. Quantum estimation is a powerful tool to infer a single parameter, as well to 
a set of parameters, up to the full reconstruction of the density matrix of an unknown 
quantum state, with or without the use of prior information [17l [18l H9] . Precision of 
any unbiased estimator is bounded by the inverse Fisher information of the probability 
distribution of the measurements outcomes, whereas the ultimate limit is written in 
term of the inverse Quantum Fisher Information (QFI). 

In quantum optical systems, homodyne measurements of field quadratures and 
Gaussian signals play a leading role. Indeed, measurement of quadratures has been 
shown to achieve phase estimation for coherent states with precision bounded by the 
(classical) Fisher information [16j . This result have been further improved by looking 
for the optimal state achieving the ultimate bound related to the QFI [20] . Among 
the pure Gaussian states, squeezed vacuum has been found to be the most sensitive 
state at fixed energy and homodyne detection |21j . Furthermore, it has been shown 
that the same signal allows optimal estimation of loss in bosonic channels [22] and 
of interaction parameters of single- and two- mode bilinear bosonic Hamiltonians [23J . 
Motivated by these results, in this paper we address optimal phase estimation by using 
Gaussian states, homodyne measurements and Bayesian estimation. We analyze the 
behavior for increasing number of measurements and show that optimality may be 
approached also with a limited number of runs upon using two-step methods acting 
on the squeezed vacuum probe and/or on the homodyne reference. Moreover, we prove 
that, in principle, the performances of double homodyne detection cannot beat the 
homodyne measurement ones, thus validating the conclusions of [21] . 

The paper is structured as follows. In Section [2] we briefly review local quantum 
estimation theory and the ultimate bounds to precision in the phase-shift estimation by 
Gaussian states. In Section [3] homodyne and double homodyne statistics are explicitly 
calculated for the phase shifted squeezed vacuum as input: this leads us to conclude 
that performances of the double homodyne detection cannot reach the limit imposed 
by QFI, while single homodyne does. Then, after describing our inference scheme, 
based on homodyne detection and Bayesian inference, the asymptotic limit for large 
number of collected data is studied in details, as well as the validity of the Gaussian 
approximation. Since the performances of this kind of inference protocol depend 
on the actual value of the (unknown) phase shift, we suggest two feasible two-step 
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adaptive methods [2H the first acting on the squeezing parameter, the other 
on the squeezing and local oscillator phases, that allow always to reach the optimal 
estimation. The results of simulated Monte Carlo experiments are reported in order 
to check convergence also for small data sample and give a quantitative meaning to 
the notion of asymptotic regime. Section [4] summarizes our results and draws some 
concluding remarks. 



2. Estimation of a phase shift 

Let us now consider a field mode undergoing a phase shift described by the unitary 
operator U(<f>) — exp(—i(/)G), with G = a^a, a and a' being the annihilation and 
creation field operators, respectively. Usually cj> itself cannot be measured and a phase 
estimation problem appears. In order to infer the value of <f> some phase-dependent 
observable X is measured and an estimator for (j>, i.e., a function of the data sample {x} 
is used. The aim of interferometry is to optimize the inference strategy by minimizing 
the uncertainty. In general, the lower bound to the variance Var[0] of any unbiased 
estimator is given by the Cramer-Rao theorem, which reads: 

Var[fl > [F^)}- 1 , (1) 

where F is the Fisher information: 

F(0)= J>(a#)[tylQgp(a#)] a , (2) 

X 

p{x\4>) being the conditional probability of obtaining the outcome x when the 
parameter has the value <f>. Since the conditional probabilities are given by p{x\4>) — 
Ti^q^Ex), Qtj) = U \4>) QqW (4>) being the quantum state of the system (actually 
depending on the initial preparation qq) and E x is the positive operator- valued 
measure (POVM) describing the measurement, Eq. ([2|) rewrites as: 

where A^, denotes the symmetric logarithmic derivative (SLD) operator: 

d<t>Q<l> = 7}(^<f>Q4> + £V>A ). (4) 

Upon using Schwartz inequality in the Hilbert space one easily shows that the Fisher 
information in Eq. ([3]) is upper bounded by the so-called quantum Fisher information 
QFI H(4>) [20], i.e.: 

F{4>) < H(4>) = Tr(^A|) . (5) 

The above equation, togehter with the Cramer-Rao theorem sets the ultimate, 
measurement-independent, bound to precision of any unbiased estimator involving 
quantum measurements. 

In order to calculate the SLD A^, we first observe that if go, and, in turn, are 
pure states, then g^ = q\ and d^g 2 ^ = {d^g^g^ + g^d^gcf,), thus, by comparison with 
Eq. ([4]), one finds A^ = 2d c f ) g l f > . More in general, we can expand go in its eigenvector 
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basis {IVVi)}) i.e., Qq = ^„Pn|V'n)(V , n| (if £>o is a pure state, then p n reduces to a 
Kronecker delta), to write: 



hk 



Then, since 



9<l>Q<t> = iy^GhkiPh ~ Pk) \^h){j>k\ (?) 

hk 

with Ghk = {iph\G\tfjk) , where \tfj n ) = U {4>)\ij) n ) , we have: 

Km | ^ = 5 (A^ |^n)(^l + . (8) 

n 

By taking the matrix elements of both sides in Eq. ((J) we obtain: 

$h\ASk) = {MA \^k) = 2iG hk ^^, (9) 

Ph + Pk 

where = U(4>)A U^((j)). As a consequence, H((f>) = Tr(^A|) = Tr(poA 2 ), i.e., the 
QFI does not depend on the value of the unknown shift cf>. The explicit evaluation of 
the QFI H = H(<t>) = H(0) leads to: 

where we used G ns — G sn . The maximum is obtained for the probe excited in a 
pure state. In this case, as described above, A^ = 2d$g g and, by substitution into 
Eq. ([5]), we obtain H = 4 AG 2 , i.e., the QFI is proportional to the fluctuations of the 
Hamiltonian G and the ultimate bound of Var[(/>] becomes: 

Var[0] = (4AG 2 )- 1 . (11) 

It is worth noticing that besides the number operator the above considerations hold 
for a general Hamiltonian generator G |23| . 

Let us now come back to the problem of estimating <fi by measurements on g^. 
Our aim is to effectively estimate the phase shift at fixed energy upon optimizing the 
measurement over detection strategies and probe states go- Of course, the ultimate 
precision is bounded by the quantum Cramer-Rao relation (jlip . which depends on 
the probe state we employ. In turn, the first stage of the optimization procedure is to 
find the best probe, which maximizes the QFI at fixed energy. We focus our attention 
onto the set of pure states and, more precisely, on Gaussian pure states, whose generic 
element is a squeezed-displaced vacuum state given by: 

Qo=D( a )S{t)\0){0\St{$D*(a), (12) 

D(a) = exp(<W — a*a) and S((,) = exp(i^a^ 2 — ^£*a 2 ), a >£ € <C, being the 
displacement and squeezing operators, respectively. In order to maximize the QFI 
we look for the state maximizing the energy fluctuations at fixed probe energy 
Tr[po^a] = sinh 2 r + |a| 2 , 

AG 2 =- sinh 2 (2r) + e 2r {Re[a\ cos <p + Im[a] sin tp} 2 

— e~ 2r {Rc[a] sintp — Im[a] cosLp} 2 , (13) 
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where we put £ = re~ 2ltp . By using Lagrange multipliers one easily finds |a| = 0: the 
maximum sensitivity is achieved when all the available energy is used to squeezed the 
vacuum. Then we have: AG 2 = isinh 2 (2r) and thus 

Var opt [0] = [2sinh 2 (2r)]-\ (14) 

which represents the ultimate bound on precision of phase-shift estimation posed by 
quantum mechanics (for Gaussian probes) [21] . Notice that Eq. (fT4|) does not depend 
on the argument ip of the complex squeezing parameter £: without lack of generality 
we will assume p — n/2. In the next Section we will show how it is possible to attain 
the ultimate precision by means of homodyne detection and Bayesian inference. 



3. Phase-shift estimation via homodyne detection and Bayesian inference 

We consider a general scheme (see Fig.[TJ where the probe state go undergoes a phase- 
shift and then the quadrature xj, is measured by homodyne detection on the outgoing 
state, OA,. 



e 





u(<t>) 




> 






probe state phase shift 



homodyne 
detection 



Figure 1. Scheme of phase estimation via homodyne detection: an input 
state go undergoes a phase shift <f>. The quadrature of the shifted state 
= U '(</>) QoU^ (<f>) is then measured by means of homodyne detection. 



The aim of our scheme is to infer the actual value 4> 01 the phase shift by processing 
the homodyne data. In order to evaluate the homodyne probability distribution we 
use the Wigner function formalism to describe our system. The Gaussian Wigner 
function associated with the state (fT2j) is (we put a — and tp = n/2): 

, . exp[— ^X T cr^. 1 X] 
Wo(X)= PL 2 J , (15) 
27r-\/Det[eroJ 

where <x = jDiag(e~ 2r , e 2r ) is the 2x2 covariance matrix. After the phase shift (see 
Fig. [J), the state is still described by a Gaussian Wigner function W$(X) of the 
form (|15p . but with covariance matrix er^ given by: 

K]n = \(e 2r cos 2 <j> + e~ 2r sin 2 0), (16) 
[<r*\ 22 = \ (e~ 2r cos 2 + e 2r sin 2 0) , (17) 
[o>]i2 = [<r$hi = ^-sinh(2r) sin(20). (18) 

At this point the quadrature = h(e measured by means of homodyne 

detection on repeated preparation of the probe state, thus obtaining a data sample 
{x}. Each outcome is distributed according to the homodyne probability distribution, 
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which can be calculated starting from the Wigner function W<f,(X) as follows: 

PtM) = ( dyW^X), (19) 



'TR 

where is a rotation matrix and X T = (x, y). Since we put tp — 7r/2, we choose to 
measure the quadrature with ip — 0. We have: 

1 ( x 2 \ 

p H (a#) = P<t>(x, 0) = exp -—J , (20) 

/27tE| v 25 V 



where: 

1 

4 

The Fisher information of the distribution (|20[) is given by 



K = - [e- 2r cos 2 </) + e 2r sin 2 0] . (21) 



F H (<f>) = Jjxp H (x\^ [d.logpuixm 2 = Sinh2( 8 2 ( ^f (20) . (22) 

Remarkably, from Eq. ([22]) we have that the Fisher information of homodyne 
distribution may be equal to the QFI upon the choice of a suitable squeezing of the 
probe state: 

r = — -log tan (23) 
or, at fixed squeezing, for a specific value of the phase shift: 

<f>n = -arcostanh 2r . (24) 

Correspondingly, the minumum variance Varjj[</>] achievable by a suitable processing 
of homodyne data mat saturate V</> to the ultimate bound (fT4"|) . 

Before going to the Bayesian inference from of homodyne data, we notice that if 
we use double-homodyne detection we have no improvement in phase-shift estimation. 
Double homodyne statistics is described by the coherent state POVM U z = 7r _1 |z)(z|, 
z e C; the probability distribution is thus given by: pT>(z\(f)) = 7r _1 |(z|J7(0)5(^)|O)| 2 
i.e.: 

exp { — \z\ 2 — tanhr Re[z 2 e 2 ^]} 

Pn{z W) = b r (25) 

7r cosn r 

where we already set £ = — r. The corresponding Fisher information reads as follows: 

F D (0) - f d 2 zp D (z\q>)[d^\og PD (z\<j))} 2 -4sinh 2 r, 

that is Fd(<P) < Fb_(4>), V(/>: the use of double homodyne detection does not bring any 
improvement of the phase-shift measurement. This result agrees with the conclusions 
of |21j , where the author considered double homodyne detection with squeezed vacuum 
as probe and an auxiliary squeezed state in the other input port. 

We stress that p(x\(f>) allows us to infer the probability of the homodyne outcome 
x once the value of 4> is assigned. In our case, the value of <fi is just the quantity 
we want to estimate and, in turn, we are interested in the conditional a posteriori 
probability distribution pu((f>\{x}) of </> given the the sample {x} = {xi, Xm} of 
homodyne data. This can be obtained by means of Bayesian inference, as we will see 
in the following. 
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3.1. Bayesian inference 

If x is the random variable associated with the outcome of the homodyne detection, 
then the Bayes' theorem states that: 

p(x\4>)p(4>) = <p{4>\x)p{x) (26) 

where p{-\-) are the conditional probabilities, p(<f>) = 2/tt is the prior assuming no 
a priori information, and p(x) the overall probability to observe x. In turn, upon 
inverting Eq. (|26p we obtain the conditional a posteriori probability p((f>\x) of 4> 
given the outcome x. After M independent homodyne measurement the a posteriori 
probability is given by 

1 M 

PM(0iw)=^n^i^)' ( 2? ) 

k—l 

AT being the normalization factor: 



N= / dct>pM{<i>\{x}). (28) 



If M > 1, then (f2Tl) rewrites as 



PM mx}) 1 Y[p(x\4) Mp W) = p{4\M) (29) 

X 

where <jf stands for the actual (unknown) value of the phase shift. In order to 
write Eq. (|29p we have used the law of large numbers and written the number of 
occurrences of the outcome x as Mp(x\(f>*). In this limit probability (|29|) can be 
explicitly calculated as follows: 

p{ct>\M) = ^exp |m Jdxp(x\r)logp(x\^ (30) 
_ 1 1 j_M^\ 

where we used log Hz — > J dx. We note that the quantity S(<fi\(f>*) = 
— ^2 x p{x\<j)*)\ogp(x\(j)) in pop may be regarded as the relative entropy between the 
two distributions [16J . In Fig. [2] the a posteriori distribution p((f>\M) is plotted for 
different values of the involved parameters as a function of (f>. It is worth noting 
that because of the asymmetric form of the distribution, a suitable estimator for the 
actual value </>* of the phase shift is given by the maximum of the distribution (its 
mode, Mode[</>]) and not to its mean: cf) = f 2 d<f) cf>p(<f)\M). This can be easily seen by 
differentiating p(<fi\M) with respect to tj>: 

dM<t>\M) = M y |M ^ W [cos(20) - cos(2r )] , (32) 
8sm(20) 



i.e., P(cj)\M) has a maximum at <fi — <fi*. However, as M increases the mode and the 
mean become the same and Eq. pip can be approximated by a Gaussian distribution 
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Figure 2. A posteriori distribution p((j>\M) for different values the number of 
data M ad squeezing parameter r. The vertical line is the actual value of the 
phase shift <j>* = 0.3. 




Figure 3. Plots of the ratio T = Sg/S^ for two values of the squeezing parameter 
r and different rj>* ■ The range of validity of Gaussian approximation strongly 
depends on the involved parameters. In particular, the less is the difference 
between </>* and the optimal phase </>h given r [see Eq. 1124 1 11 . the larger is the 
range of validity of this approximation (in the plots we have r = 0.2 — > (f'Yi = 0.59 
and r = 0.6 -> H = 0.29). 



with mean <\>* and variance Y? g given by: 



T 2 - 



1 



d 2 p(<j)\M) 



1 

M 



p{4>*\M) d(j> 2 

1 d 2 p{x\ 



E 



p(x\(j)*) d<j) 2 



MF(<f>*y 



(33) 
(34) 



where we substituted Eq. (|3"D|) into Eq. (|3l?l) and F '{<$*) is the Fisher information of the 
probaility distribution p(x\4>*). The factor M^ 1 follows from taking the data sample 
as a collection of M mutually independent measurements, which, indeed, leads to an 
ensemble average over M different copies of the system. Finally, we notice that, as 
one may expect, the variance and, thus, the precision of the estimation depends on 
the true value 4>* itself. 

Overall, the Bayes estimator is asymptotically unbiased and efficient, i.e., the variance 
Var[0] saturate the Cramer-Rao bound of Eq. this is a consequence of the 

asymptotic normality of the a posteriori distribution (Laplace-Bernstein-von Mises 
theorem) [27} [28] . However, two questions arises. The first concerns the range of 
validity of the Gaussian approximation, which depends on both cj>* and the squeezing 
parameter r. This aspect is illustrated in Fig.[3]where we plot the ratio T = Sg/S^, E B 
being the variance of the asymptotic distribution p{<f>\M). For fixed r, one finds that 
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Figure 4. Plot of the ratio R(r) = M£g(r)/Var op t[0*] (see text) as a function 
of r and for <f>* = 0.3. The vertical dashed line indicates r op t- 

the less is the difference between (f>* and the optimal phase <pH given r [see Eq. (]24[)]. 
the larger is the range of validity of this approximation. On this observation is also 
based the two-step adaptive method we will describe below. The second question is 
whether the Bayes estimator may saturate also the quantum Cramer- Rao bound, i.e., 
whether the Fisher information of p(x\4>) may be equal to the QFI, thus leading to 
phase-shift estimation with precision at the ultimate quantum limit. As what concerns 
this point we notice that, being the variance of Bayes estimator dependent on the true 
value of the phase shift, some kind of feedback should be unavoidably involved. In the 
following we will describe two possible adaptive mechanisms, acting on the squeezing 
parameter of the probe or on the homodyne local oscillator and squeezing phase, 
respectively. 

3. 2. Examples of two-step methods to achieve ultimate precision 

Adaptive methods for Bayesian estimation allow to always attain the ultimate bound 
on precision and have been investigated in the case of large ensembles and qubit 
systems 130]. Here we propose two realistic and feasible setups exploiting the 
interferometric features of homodyne detection. 

The first scheme is based on the fact that the variance S^(r) may achieve the 
optimal value Af~ 1 Var opt [(f)*] of Eq. (fl4f employing a squeezed vacuum probe with 
parameter r opt = — \ log tan (f>* . Of course, setting r = r opt requires the knowledge 
of the actual (unknown) value of the phase shift. However, one may obtain a rough 
estimate of </>* upon building the distribution p(<f)\M') with a fraction of the M, taking 
its maximum (Mode[0]) and then modify the probe state, tuning its squeezing to r opt . 
In Fig. H we show the ratio R(r) = MS^(r) /Var op t[0*] for the case 4>* = 0.3: the 
smooth behavior of R(r) ensures the convergence of the above mechanism. Tuning 
the squeezing parameter, however, could be a challenging task. On the other hand, also 
when r and, thus, the energy are fixed, it is possible to achieve the optimal variance 
by tuning the squeezing phase (p of the probe state or the phase ip of the homodyne 
quadrature. In fact, previously we set ip = ir/2 and ip = 0; if, on the contrary, we 
assign to these phases the generic values ip and ip, then we should simply apply the 
following change of variable in all the previous equations: cj> — * <f>+ Up — ip — ^J, that 
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Figure 5. Bayesian estimation of the phase-shift from Monte Carlo simulated 
homodyne measurements: on the left the ratio A = <j>/<t>* an d on the right 
V = \/Var [4>]/ Var op t [<t>] (right). The solid lines denote results obtained with the 
adaptive method acting on the homodyne local oscillator and squeezing phases; 
the dashed lines are obtained without the adaptive method. We set r = 0.6 and 
4>* = 0.7. In the case of the adaptive method, we used [_3a/M] of the M data to 
argue the phase-shift rough estimate (dotted line), then the left homodyne data 
are processed to assess and Var [<j>]. In both the experiments we use the same 
total number M of data. Lower panel: the same as in the top panel for r = 0.3. 



is a translation of <f> by the amount <p — if) — ^ . Since the optimal angle </>h at fixed r 
is given by Eq. (|24p , optimality is always achieved by choosing: 

7T 

(p-ip = ^H-^+-. (35) 

As described above, we may obtain a rough estimate of <ff by taking the maximum 
(Mode[0]) of p((f)\M), and, thus, we can tune the quantity ip — t/j by means of adaptive 
control on the homodyne detection and/or the probe state, whose energy does not 
depends on ip. 

In order to confirm convergence also for small data sample, we performed a set 
of Monte Carlo simulated experiments with the latter adaptive scheme. The results 
are shown in Fig. [S]for r = 0.6 and 0* = 0.7 (upper panel) and with reduced energy 
r = 0.3 (lower panel). In the experiment without adaptive method the whole sample 
of M homodyne data, obtained as described in the first part of this section, is used 
to estimate 4> and Var[0] (dashed lines in Fig. [5|). With the adaptive scheme (solid 
lines), iV r = [3VMJ of the M data sample are used to argue the phase-shift rough 
estimate, then the phase difference ip — ip is tuned according to Eq. (|35p and the left 
homodyne data are processed to assess <j> and Var [(/>]. Each point in Fig. \5\ corresponds 
to the average over 20 repetitions. Of course, the effectiveness of the adaptive method 
depends on the value of the rough estimate: in this view, an increasing number of the 
outcomes devoted to the rough estimation, as the data sample becomes larger, allows 
the reduction of the Var[</>] fluctuations, as one may verify, for example, by using a 
fixed value for N T . It is worth to note that in our simulations the rough estimate is 
obtained as Mode[</>], whereas the mean (j> is used for the final results: this is justified 
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for the small N T considered in the rough estimate and the larger number of the final 
estimation (the error introduced by this choice does not sensitively affect our results, 
as we verified also assessing the Pearson skewness coefficient \<f> — Mo de [(f)] | / v/Var [(f)]). 

4. Conclusions 

In this paper we have shown how Bayesian inference techniques represent useful tools 
for phase estimation. Our analysis is based on homodyne detection with squeezed 
vacuum as a probe state, and Bayesian post-processing to infer the phase shift. In 
the asymptotic limit of a large number of measurements, our scheme saturates the 
Cramer- Rao bound to precision, i.e., the variance of the phase shift achieves the 
lower bound imposed by the inverse Fisher information. Moreover, we have shown 
that optimality may be approached also with a limited number of measurements 
by means of two-step methods acting on the squeezed vacuum probe and/or on 
the homodyne reference. These have been investigated by means of Monte Carlo 
simulated experiments, which show excellent results also in the case of small data 
samples. Our results, together with the recent advances in homodyne detection |31j 
lead us to conclude that the estimation protocol described in our paper may be suitable 
for experimental investigation, opening the way to information technology based on 
Gaussian states and phase encoding. 
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